#ifndef NAKAGAMI_PARAM_ESTIMATION_H_INCLUDED
#define NAKAGAMI_PARAM_ESTIMATION_H_INCLUDED


 const static double bvalues[] = {	/* Bernoulli Numbers */
	 1.00000000000000000e+00,
	 1.66666666666666667e-01,
	-3.33333333333333333e-02,
	 2.38095238095238095e-02,
	-3.33333333333333333e-02,
	 7.57575757575757576e-02,
	-2.53113553113553114e-01,
	 1.16666666666666667e+00,
	-7.09215686274509804e+00,
	 5.49711779448621554e+01,
	-5.29124242424242424e+02,
	 6.19212318840579710e+03,
	-8.65802531135531136e+04,
	 1.42551716666666667e+06,
	-2.72982310678160920e+07,
	 6.01580873900642368e+08,
	-1.51163157670921569e+10,
	 4.29614643061166667e+11,
	-1.37116552050883328e+13,
	 4.88332318973593167e+14,
	-1.92965793419400681e+16
    };
const static int K=14;
const static double x0=7.00000000000000000;
const static double Rx=-2.00000000000000000e-2;
float pow_func(float x, float y)
{

	float temp1, temp2;
	 temp1=(float)x/fabs(x);
	 temp2=y-floor(y);
	 if(temp1<0 && temp2>=0)
		{	
			return (-1)*pow(fabs(x), y);
			

		}	
	else
		return pow(x,y);	



}
int getInt(double num)
{
	double temp;
	
	temp=(num-floor(num));
	if(temp<0.5f)
		return (int)floor(num);
	else
		return (int)floor(num+1);
}



double pshifunction(int n, float z)
{
		
	int m,  r_indx;
	double t1, interm_val1, interm_val2,interm_val3, f_val, intpart1,intpart2, temp1; 
       
	if(x0>0 && x0<=z)
	{
		 if(n==0)
		{
			t1=-log(z);
		}
		else 
		{
			t1=(double)factorial((double)(n-1))/pow(z,n);
			
		}
			
		interm_val1=(double)factorial((double)n)/(2*pow(z, (n+1)));
		interm_val2=0.0;
		for(r_indx=1;r_indx<=K;r_indx++)
		{
				
	   interm_val2=interm_val2+(bvalues[r_indx]*(((double)factorial((double)(2*r_indx+n-1)))/((factorial((double)(2*r_indx)))*pow(z,(2*r_indx+n)))));
			//printf("%0.30f\n", interm_val2);				
		}
							
				f_val=pow(-1,n-1)*(t1+interm_val1+interm_val2);
			return f_val;
		}
		else if(z>=0 && z<x0)
		{
			temp1=modf(x0, &intpart1);
			temp1=modf(z, &intpart2);
			m=(int)(intpart1-intpart2);
				
			interm_val3=pshifunction(n, z+m);
			interm_val2=0.0;
			for(r_indx=0;r_indx<=m-1;r_indx++)
			{
				interm_val2=interm_val2+(double)1/pow((z+r_indx),(n+1));
				
			}
			f_val=interm_val3-(pow(-1,n)*factorial((double)n)*interm_val2);
			return f_val;

		}

}
float calc_mean(float **img, int height, int width)
{
	int row, col;
	float mean=0.0;
	for(row=0; row<height; row++)
	{
		for(col=0;col<width; col++)
		{
			mean=mean+img[row][col];
		}
	}
	mean=(float)mean/(height*width);
	return mean;

}
float calc_std_dev(float **img, float  mean, int po, int height, int width)
{
	int row, col;
	float var=0.0, **dev;
	dev=alloc2Dfloat(height, width);
	for(row=0; row<height; row++)
	{
		for(col=0;col<width; col++)
		{
						
				dev[row][col]=pow_func((img[row][col]-mean), po);
				
			
		}
	}
	for(row=0; row<height; row++)
	{
		for(col=0;col<width; col++)
		{
			var=var+dev[row][col];
		}
	}
	var=(float)var/(height*width);
	return var;
}
float newton_method(float initial_guess, float gx_val1, float gx_val2)
{
	int iter=0; 
	float x,  f_x, f_dash, psi_1, psi_2, psi_3, error, f_x_new,prev_x;
	x=initial_guess;
		psi_1= pshifunction(1, x);
		psi_2= pshifunction(2, x);	
		psi_3= pshifunction(3, x);
		f_x_new=psi_2*pow_func(gx_val1,1.5)-pow_func(psi_1,1.5)*gx_val2;
	do
	{
        prev_x=x;
	   	f_x=f_x_new;
		f_dash=psi_3*pow_func(gx_val1,1.5)-1.5*psi_2*gx_val2*pow_func(psi_1,1.5);
		//printf("\n%f\t%f\t%f\t%f\t%f\t%f\t%f\t%f\t%f\n",psi_1,psi_2,psi_3, f_x, f_dash, gx_val1, gx_val2,x);
		x= x+(float)f_x/(f_dash);
		
		psi_1= pshifunction(1, x);
		psi_2= pshifunction(2, x);	
		psi_3= pshifunction(3, x);
		f_x_new=psi_2*pow_func(gx_val1,1.5)-psi_1*gx_val2;
		error=sqrt((f_x_new-f_x)*(f_x_new-f_x));
		iter++;
		//printf("%d\t%f\t%f\n",++iter,x, error);
	}while(error>0.05f && iter<5 );	
  
	return x;

}

float m_estimation(float **img, int height, int width)
{
  float mean_im, std_dev1, fac1, fac2, estimated_m;
  int  c_indx;
  mean_im=calc_mean(img, height, width);
  std_dev1=calc_std_dev(img, mean_im,2, height, width);
  fac2=calc_std_dev(img, mean_im, 3,height, width);
  fac1=pow_func(std_dev1, 1.5);
	printf("fac1:%f\tfac2:%f\n",fac1, fac2);
  estimated_m=newton_method(1, fac1,fac2);
  return estimated_m>0.5 ? estimated_m: 0.5;
}

float s_estimation(float **img, float m, int height, int width)
{
	float mean, std_dev1, std_dev2, psi_1, psi_2, estimated_s, Num, Den;
	mean=calc_mean(img,height, width);
	std_dev1=calc_std_dev(img, mean, 2, height, width);
	std_dev2=calc_std_dev(img, mean, 3,height, width);
	psi_1=pshifunction(2,m);
	psi_2=pshifunction(1,m);
	Num=fabs(std_dev1*psi_1);
	Den=fabs(std_dev2*psi_2);
	
	estimated_s=0.5*fabs((float)Num/(Den+1));
	return estimated_s;
	
}

float omega_estimation(float **img, float m, float s, int height, int width)
{
		float mean, psi1, temp1, estimated_omega;
		mean=calc_mean(img, height, width);
		psi1=pshifunction(0, m);
		temp1=exp(mean-(float)psi1/(2*s+1));
		
		estimated_omega=m*pow_func(temp1, (2*s));
		
		return estimated_omega;
}	
float m_estimate1D(float *vec, int length)
{
	int indx;
	float mean=0, var1=0, var2=0, m;
	for(indx=0;indx<length;indx++)
	{
		//printf("vector: %d\t%f\t %f\n", indx, vec[indx], log(vec[indx]));
		mean=mean+log(vec[indx]+1);
	}
	mean=(float)mean/length;
	for(indx=0;indx<length;indx++)
	{
		var1=var1+pow_func((log(vec[indx]+1)-mean),2);
		var2=var2+pow_func((log(vec[indx]+1)-mean),3);
	}
	var1=(float)var1/length;
	var2=(float)var2/length;
	var1=pow_func(var1,1.5);
	
	m=newton_method(1, fabs(var1),fabs(var2));
	m=m>0.5 ? m: 0.5;
      	return m;

}

float s_estimate1D(float *vec, float m, int length)
{
	int indx;	
	float mean=0, var1=0, var2=0, psi_1, psi_2, estimated_s, Num, Den;
	mean=0;
	for(indx=0;indx<length;indx++)
	{
		mean=mean+log(vec[indx]+1);
	}
	mean=(float)mean/length;
	for(indx=0;indx<length;indx++)
	{
		var1=var1+pow_func((log(vec[indx]+1)-mean),2);
		var2=var2+pow_func((log(vec[indx]+1)-mean),3);
	}
	var1=(float)var1/length;
	var2=(float)var2/length;
	psi_1=pshifunction(2,m);
	psi_2=pshifunction(1,m);
	Num=fabs(var1*psi_1);
	Den=fabs(var2*psi_2);
	//printf("Den:%f\tstd_dev2:%f\t\n", Den,var2);
		
	estimated_s=0.5*fabs((float)Num/(Den+0.0001));
		
	estimated_s=estimated_s>0.5 ?estimated_s: 0.5;
	return estimated_s;
	
}

float omega_estimate1D(float *vec, float m, float s, int length)
{
		int indx;		
		float mean, psi1, temp1, estimated_omega;
		mean=0;
		for(indx=0;indx<length;indx++)
		{
			mean=mean+log(vec[indx]+1);
			
		}
		mean=(float)mean/length;
		psi1=pshifunction(0, m);
		temp1=exp(mean-(float)psi1/(2*s));
		
		estimated_omega=m*pow_func(temp1, (2*s));
		
		return estimated_omega;
}	


#else
#endif


/*int main(int argc, char *argv[])
{
	double num, f_val, res;
	num=atof(argv[1]);
	f_val= pshifunction(2, num);
	res=newton_method(num, atof(argv[2]));
	//printf("The factorial value is %f",factorial(2));
	printf("The psi value is %f",f_val);
	return 0;
}*/
